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Abstract 

Recently, optimization techniques are increasingly being used for per- 
forming nonlinear structural analysis. The development of element-by-element 
(EBE) preconditioned conjugate gradient (CG) techniques is expected to extend 
this trend to linear analysis. Under these circumstances the structural 
design problem can be viewed as a nested optimization problem. The present 
paper suggests that there are computational benefits to treating this nested 
problem as a large single optimization problem. That is, the response 
variables (such as displacements) and the structural parameters are all 
treated as design variables in a unified formulation which performs simul- 
taneously the design and analysis. Two examples are used for demonstration. 

A seventy- two bar truss is optimized subject to linear stress constraints 
and a wing-box structure is optimized subject to nonlinear collapse constraints. 
Both examples show substantial computational savings with the unified approach 
as compared to the traditional nested approach. 

Introduction 

Structural optimization was initially based on the calculus of 
variations. A typical problem was solved by obtaining the Euler-Lagrange 
optimality differential equations and solving them simultaneously with the 
differential equations of the structural response. This approach is still 



used today for the optimization of individual structural elements such as 
beam-columns (Ref. 1). However, for built-up structures modeled by finite 
elements, a nested approach is typical. Resizing rules based on optimality 
criteria require that the structural response be calculated repeatedly for 
each set of trial structural design variables (see, for example. Ref. 2). 
This preference for the nested rather than the simultaneous approach is 
probably due to the simplicity of the structural resizing rules which are 
possible when the structural response is known. This simplicity contrasts 
with the difficulty of solving the large systems of nonlinear algebraic 
equations which are obtained from a simultaneous formulation. 

In the last twenty-five years direct search methods have been gaining 
ground as the standard for structural optimization. These techniques are 
commonly used in a nested approach. That is, the structural analysis 
equations are repeatedly solved during each design iteration. Part of the 
reason for the popularity of the nested approach is that the structural 
analysis equations are solved by techniques which are quite different than 
those used for the design optimization. An exception is the design of a 
structure subject to constraints on its collapse load. There, the analysis 
problem ("limit analysis") is often approximated as a linear program and 
solved by the simplex method. The structural design problem in that case 
("limit design") is easily formulated as a single linear program with the 
element forces and structural parameters both treated as design variables 
(Ref. 3, for example). 

In the late sixties, Schmit, Fox and their coworkers (Refs. 4-6) tried 
to unify the treatment of structural analysis and design by employing 
conjugate gradient (CG) minimization techniques for solving linear structural 
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analysis problems. They found that the optimization methods were not 
competitive with the traditional direct Gaussian elimination techniques. 

More recently techniques for unconstrained minimization have become 
more efficient and their application to structural analysis has become 
more feasible (e.g. ; Ref. 7). The emergence of the preconditioned 
conjugate gradient techniques (e.g., Ref. 8) and the element- by-element 
(EBE) formulations of Hughes and coworkers (Ref. 9) make CG methods 
particularly attractive for structural analysis. 

In view of the increasing use of optimization methods for structural 
analysis, there is merit in considering again the simultaneous approach to 
analysis and design. When optimization techniques are used for structural 
analysis the design problem becomes a nested optimization problem. Such 
nested problems can, indeed, often benefit computationally from a single 
level treatment that disregards the nested structure (Ref. 10). The present 
paper investigates the use of the simultaneous analysis and design approach 
for linear and nonlinear problems. An EBE preconditioned CG method is 
applied to a linear analysis problem, and Newton's method is used for design 
subject to a nonlinear collapse load constraint. 
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Linear Structural Analysis and the Conjugate 
Gradient Method 

In the linear case the finite element discretization of a structure 
typically leads to a system of linear equations 

R s KU - F = 0 (1) 

where K is a stiffness matrix, U a displacement vector and F a load vector. 
Eq. (1) is usually solved by Gaussian elimination, but it can also be 
solved by minimizing the function 


f(U) = 1/2 (U-U) T K(U-U) (2) 

where U is the solution of Eq. (1). Clearly the minimum of f is obtained 
for U = U. The function f(0) may be minimized even though U is unknown because 
its gradient G may be calculated without using U 

G = vf = K(O-U) = KU-F (3) 

Applying the conjugate gradient (CG) method to the minimization of f is 
often very expensive because the second derivative matrix of fj which is 
usually has a high condition number. This problem can be remedied if we 
have an approximation B to K. Then B can be used to "precondition" the 
problem. The preconditioned CG algorithm is summarized in the Appendix. 

While any reasonable approximation to the inverse of K would solve the 
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ill conditioning problem, the method is economical only if B" 1 is cheap to 
calculate compared to K \ When K is very poorly banded one approach is 
to obtain B by an incomplete factorization of K (e.g., Ref. 11). Another 
approach which is attractive is the element-by-element (EBE) technique of 
Hughes (Ref. 9). Reference 9 presents several formulations for obtaining B, 
and the one used herein is given in the Appendix. 

Simultaneous Analysis and Design: Linear Case 

The structural design problem may be formulated as 
find X to minimize m(X) 

subject to gj(U,X) L 0 j = l,...m (4) 

where m is an objective function, X a vector of design parameters and g. 

J 

constraint functions such as stress and displacement constraints. In addition, 
the displacement vector U must satisfy Eq. (1) where both K and F may depend 
on X. The optimization problem is usually solved by repeatedly calculating 

U and its derivatives with respect to the components of X, x.. U is calculated 

dU ^ 

from Eq. (1) and -j^- is calculated either by differentiating Eq. (1) or by 

J 

finite differences. Based on U and its derivatives the constraint functions 
9j and their derivatives can be evaluated and a numerical optimization technique 
can use this information to improve X. 

This approach is wasteful because the displacement vector U is calculated 
exactly for each trial design X. A simultaneous approach treats X and U equally 
as design variables and solves the following expanded problem 

find X,U to minimize m(X) 

Subject to 

gjtt.tO^O j = l,...,m 
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and 


R = KU - F = 0 (5) 

The equations of equilibrium are treated here as equality constraints. 

The optimization technique used here to solve problem (5) is a penalty 

function technique. An application of a standard penalty technique replaces 
(5) by 


minimize <f>(X,U,r) . = cm(X) + r z p[g.(X,U)] + ~ R T R 


m _ . c i t 

J n » r 


for 


j=l 


/ r 


r r^ , rg , . . . 


where 


where 

p(g) = ^-[(g/g 0 ) 2 - 3g/g Q + 3] ( 7 ) 

is an extended interior penalty function (Ref. 12) with g Q being a 
transition parameter. The constants c and c-j are chosen to balance the 
contribution of the objective function, inequality constraints and equality 
constraints to <j> (see Ref. 13 for details). If the minimization of <p is 
accomplished by the CG method then the R^R term presents a problem 

R T R = (KU-F) T (KU-F) (8) 

so that the second derivative matrix of R T R with respect to U is K T K. 
Unfortunately, the condition number of K T K is the square of that of K so 
that the optimization of $ by the CG method would proceed extremely slowly. 
It is tempting to replace the term R T R by the function f defined by Eq. (2). 
However, while the derivatives of f with respect to the displacement can be 
calculated (see Eq. (3)) without knowing the exact displacement vector U, 
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the derivatives of f with respect to X cannot. 

It is due to the ill-conditioning of <p that the idea of simultaneous 
analysis and design has not been pursued in the past. The solution proposed 
here is to replace the term R T R by R T B _1 R so that 

*(X.U,r) = cm(X) + r? p[g.(X,U)] + ^k T B -1 R (9) 

j=l J yfr v ' 

Because B is an approximation to K" 1 the R T BR term has similar conditioning 

as f. 

Another way of avoiding the ill-conditioning due to the R T R term is to 
use Newton s method instead of the CG method. This option is explored in the 
next section. 
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Simultaneous Analysis and Design: Nonlinear Case 

The formulation described in the previous section can also be used in 
the case where R is a nonlinear function of U. In that case B would be an 
approximation to the Jacobian of R. To broaden the scope of the dis- 
cussion, a different nonlinear formulation based on collapse techniques 
(Ref. 13) is presented below. This formulation contains a combination of 
linear and nonlinear analyses and demonstrates the option of solving the 
nonlinear analysis problem simultaneously with the design problem while treat 
ing the linear analysis sequentially. 

Reference 13 considers a structure which is subject to some load 
conditions when it is intact, and to other load conditions when it is damaged 
For the sake of simplicity only one load case in each condition is treated 
here. In its undamaged condition the structure is subject to displacement 
constraints 


stress constraints 

g di ( x » u ) 1 0 

i - 1 , . . . , n^ 

(10) 

and buckling constraints 

g si (x,u) > o 

i - 1,..., n s 

(11) 

The displacement U can be 

g bi (x,u) i o 

calculated from the 

i = 1 , . . . , n b 
linear analysis 

(12) 


KU = Fj 


(13) 


where Fj is the load applied to the undamaged structure. For the damaged 
structure, large deformation and post-buckling response can be tolerated as 
long as the structure does not collapse. In this case, an approximate 
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calculation of the collapse load can be performed by assuming that the 
equations of compatibility can be disregarded. That is, we assume that 
after yielding or buckling the internal loads begin to redistribute, with 
yielded and buckled elements becoming "soft" and undergoing large deformations. 
The situation is idealized here by assuming zero post-buckling or post-yielding 
stiffness. That is, the yielded or buckled element continues to carry the load 
that it carried at the onset of buckling or yielding, and that load does not 
increase with additional deformation. The collapse load is reached when no 
amount of internal load redistribution can balance the applied loads. To take 
advantage of the above assumption for calculating the collapse load, the 
element forces are used as the unknowns in the equations of equilibrium instead 
of the displacements. That is, the equation of equilibrium are written as 

ET = fF D (14) 

where E is a matrix which depends only on the initial geometry of the structure, 
T is a vector of element internal loads, F p is the design load for the damaged 
structure, and f is a safety factor. For statically indeterminate structures 
the matrix E is rectangular and it is not possible to determine T uniquely 
from Eq. (14). Instead, based on the assumptions that we made, the structure 
will not collapse if it is possible to find a solution to Eq. (14) for f = 1 
such that the stress and buckling constraints are satisfied. These constraints 
are rewritten in terms of T as 

9 si ( T )-0 i = l,...,n s (15) 

9 bi (T)iO i = 1 n b (16) 

The constraint that the structure does not collapse at the applied load 
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may be written as 


*max- 07) 

where f max is the maximum value of f that may be obtained by selecting T 
that satisfies Eqs. (14-16). 

In Reference 13 the design problem was posed as 
find X to minimize m(X) 
subject to 


g di (u) > o 

i = 1 , . . . , n^ 


g si (u) >o 

i - 1,..., n g 


g bi (u) > 0 

i = 1 .... .n^ 

(18) 


and 


f > 1 

max - 


where U was obtained by a direct solution (Gaussian elimination) of Eq. 

(13) and f max was obtained by solving a maximization problem. The calcu- 
lation of the displacement vector U and the collapse margin of safety 
f max re P resents the structural analysis and had to be repeated many times 
as the optimization with respect to X proceeded. However, the calculation 
of f max was by far more cost 1y than the calculation of U. For this reason, 
the technique of simultaneous analysis and design is applied here only to 
the collapse constraint. That is, the design problem is rewritten as 


find X,T to minimize 

m(X) 


subject to 




g di (u) >o 

i = 1 


g sj (u) iO 

i - 1 


g b 1 (u) iO 

i = 1 


ET = fF r 
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g si (T) > 0 i = l,...»n s 

9 bi (T) lO i - 1 , . . . , n^ 

f>l (19) 

where U is calculated from Eq. (13). The optimization problem (19) is 
solved by a penalty function approach by minimizing an augmented objective 
function 

n d " s 

<t>(X,T,r) = cm(X) + r p[g di (U)] + r ? =1 {p[g $i (U)] + p[g si (T)]} 
n b c 

+ r LJP [ 9b1 (U)] + PE9bi (T)] l + /F (ET-fF D ) T (ET-fF D ) + rp(f-l) 

( 20 ) 

To overcome the ill-conditioning coming from the penalty due to the equations 
of equilibrium, Newton's method rather than a preconditioned CG method was 
employed. This was motivated by the fact that unlike the matrix K, the 
matrix E is not a function of the design variables so that second derivatives 
are easier to calculate. The second derivatives of the penalty terms for 
the undamaged structure were approximated by using only first derivatives of 
the constraints (see Ref. 12). 


11 



RESULTS AND DISCUSSION 
Linear Analysis 

The example used for demonstrating the usefulness of an EBE pre- 
conditioned CG method for the linear case is a 72 bar truss shown in 
Figure 1. The loading and the stress allowables are given in Table 1. 

The difference between the regular CG method and the EBE preconditioned 
CG method is shown in Figure 2 which shows the convergence of the stress 
in member 1. It is clear from Figure 2 that the convergence of the pre- 
conditioned CG method is much faster than of the regular CG method. 

Next, the truss was optimized subject to stress constraints. The ' 
cross-sectional areas of the 72 members were the design variables. The 
optimization was performed three different ways. First the traditional 
sequential approach was employed using a standard truss finite element 
analysis to calculate stresses and the NEWSUMT optimization program (Ref. 14) 
to optimize the cross-sectional areas. The derivatives of the constraints 
were calculated by finite differences. The second optimization was carried 
out by a CG package based on Beale's restarted CG algorithm (Ref. 15) 
applied to the penalty function formulation of Eq. (6). The design vari- 
ables were both the 72 cross-sectional area and the 48 nonzero displace- 
ment components. Finally the same optimization was repeated with the 
modified formulation of Eq. (9) where the EBE approximate matrix B was 
calculated for the initial design and not updated as the design changed. 

The results of the three optimizations are summarized in Table 2. As can 
be seen from Table 2 the use of the EBE-generated approximate inverse 
reduced the number of constraint evaluations and CPU time by an order of 



magnitude. The comparison with the traditional sequential optimization 
approach is more difficult because its cost is inflated by the use of 
finite difference derivatives, and the optimization method is different. 
Even so, it does appear that the simultaneous optimization approach could 
be expected to be better than the traditional sequential approach. 

The results of the 72 bar truss were obtained for a single load case. 
When multiple load cases need to be considered the CG method is at a 
disadvantage compared to direct methods. In such cases the simultaneous 
approach could rely on Newton's method. The use of Newton's method is 
demonstrated now for a nonlinear problem. 

Nonlinear Analysis 

Reference 13 discusses the design for damage tolerance of the wing- 
box structure shown in Figure 3. The box is 3.56 m (140 in.) long, 2.24 m 
(88 in.) wide, and 38 cm (15 in.) deep. As shown in Fig. 3, the wing box 
is clamped at the root and a variable load is applied at the tip. The 
loads applied at the four tip nodes are 8163, 17,000, 17,000, and 34,000 
kg (18,000, 37,500, and 75,000 lb). The upper and lower wing-skin panels 
are modeled by membrane elements, and the webs of the ribs and spars are 
modeled by shear web elements. The spar caps and vertical posts at the 
rib- spar intersections are modeled by rod elements. The wing is assumed 
to be made of 7075 aluminum alloy. 

The finite element model employs a single quadrilateral membrane 
element to represent a skin panel bounded by ribs and spars. The model 
has 32 grid points and 75 finite elements. Design variables are distri- 
buted so that there is one design variable for the thickness of each cover 
skin between adjacent ribs (i.e., constant chordwise distribution), one 
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design variable assigned to the thickness of each spar web between ribs, 
one design variable assigned to the area of each spar cap between ribs, 
and one design variable assigned to the thickness of each rib. The cross- 
sectional areas of the vertical posts at the rib-spar intersections are 
held constant. The total number of design variables is 45. The design 
constraints are stress constraints [503 MPa (73,000 psi) allowable stress, 
using the von Mises yield criterion], buckling constraints for panels 
and shear webs, and side constraints. The side constraints were 0.5 mm 
(0.02 in.) minimum gage for skin panels, 0.65 cm 2 (0.1 in. 2 ) minimum area 
for spar caps, and 3.9 cm 2 (0.6 in. 2 ) maximum area for spar caps. 

The computation time with the traditional sequential approach was 
7000 CPU sec. on a CDC Cyber-173 computer. The simultaneous approach 
required only 1600 CPU sec. The final designs had similar design variable 
distribution with less than 1% difference in mass. 

CONCLUDING REMARKS 

Nonlinear structural analysis is usually performed by iterative 
methods. Recently iterative algorithms are becoming competitive even for 
linear analysis. The present paper investigated the possibility of inter- 
facing the analysis iterations with design optimization iterations and 
combining the two in a single optimization problem. 

For linear problems a seventy- two- bar truss example was used to show 
that the Element-by-Element approximate inverse of the stiffness matrix can 
be used to speed the convergence of the simultaneous analysis-design con- 
vergence by an order of magnitude. 

For a nonlinear wing-box damage- tolerance example employing collapse 
techniques, the simultaneous approach reduced the computation time by 
better than a factor of four. 
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Appendix - The EBE Preconditioned CG Method 


Step 

Step 

Step 

Step 

Step 


Step 

Step 

Step 


The 

flow chart for solving the 

linear system KU = F is 

1: 

Initialization m = 0, U 0 = 

0, R„ = F, P = Z = B -1 R 

OOO 0 

2: 

u = R m T Z /P m T KP m 
m m m m m 


3: 

^ntf-l ^m + “m P m 


4: 

R m+1 = R m ' a m KP m 


5: 

Convergence check 



HVill ? 



Yes: Return 



No: Continue 


6: 

Z »*l " 


7: 

B m " *Wl Z m*l Z ^m Z m 


8: 

P ntf-1 " Z m+1 + & m P m 



Increment m and go to step 2 

The matrix Bis an approximation of K obtained (Ref. 9) as 

B = W 1/2 C W 1/2 


where W is a diagonal scaling matrix taken here to be the diagonal of the 
matrix K. The matrix C is a product of factored element matrices given 
below. For each finite element we generate the factorization 



where 


K e = W _1/2 (K e - D e )W" 1/2 

0 

K is the element stiffness matrix and D e is a diagonal matrix composed of 
the diagonal terms of K . The matrix C is then given as 



r n el 1 


r n el 1 


‘1 

c = 

» L P 


n D p 


I' (Lp) T 


, e=l 


_e=l . 


- e=n el 


and n e -j is the number of finite elements. The factorization of C is 
performed completely at the element level. 




Table 1 

Definition of 
72-Bar Space Truss 
(U. S. Customary Units) 

T 

Material 

: Aluminum 

a 

Young's Modulus 

: E 

= 10 7 psi 


Specific mass 

• p 

= 0.1 lbm/in 3 


Allowable stress 

°a 

= + 25,000 psi 


Minimum area 

: 0 < L > 

= 0.1 in 2 


Uniform initial area 
Nodal loading 

: 

= 1.0 in 2 



Load components (Ibf) 


5,000 


5,000 


-5,000 



-5,000 

-5,000 

-5,000 











Table 2: Comparison of Sequential and Simultaneous 

Analysis-Design Optimization for 72 Bar Truss 


Approach 

Solution 

Method 

No. of 
r Values 

No. of 

Constraint 

Evaluations 

CPU Time 
IBM 3081 
(sec) 

Final 

Weight 

(lb) 

Sequential 

Newton 

14 

4733 (4464 
for derivatives) 

185.9 

95.6 

Simultaneous 

CG 

12 

13347 

150.0 

96.7 

Simultaneous 

EBE-CG 

12 

665 

8.8 

96.7 
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Figure 2: Convergence of Stress in Member 1 of 72-Bar Truss 





DAMAGE CONDITIONS 



Figure 3: Wing Box Structure 
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